setwd('/data1/bsi/bioinf_int/s106381.borawork/data_HTHGU/converted_plus2_customcdf_frma/convertedhgu_replicates')


library(preprocessCore)
library(MASS)
library(affy)
library(simpleaffy)
library(frma)
data <- ReadAffy()


data@cdfName<-"newone1hgu133ageneccds"

frma1 <- function(object, background="rma", normalize="quantile", summarize="robust_weighted_average", input.vecs=list(normVec=NULL, probeVec=NULL, probeVarBetween=NULL, probeVarWithin=NULL, probesetSD=NULL), output.param=NULL, verbose=FALSE){

  if(!is(object, "AffyBatch")) stop(paste("argument is", class(object), "frma requires AffyBatch"))

  if(!background %in% c("none", "rma")) stop("background must be either none or rma")
  if(!normalize %in% c("none", "quantile")) stop("normalize must be either none or quantile")
  if(!summarize %in% c("median_polish", "average", "median", "weighted_average", "robust_weighted_average", "batch")) stop("summarize must be one of: median_polish, average, median, weighted_average, robust_weighted_average, batch")

  cdfname <- cleancdfname(cdfName(object))
  platform <- gsub("cdf","",cdfname)

  if(summarize == "median_polish") output <- frmaMedPol(object, background, normalize, input.vecs, verbose)
  if(summarize %in% c("average", "median", "weighted_average", "robust_weighted_average")) output <- frmaRobReg1(object, background, normalize, summarize, input.vecs, output.param, verbose)
  if(summarize == "batch") output <- frmaBatch1(object, background, normalize, input.vecs, output.param, verbose)

  if(is.null(output.param)){
    e <- new("ExpressionSet", exprs=output$exprs, annotation=annotation(object))
  } else if("stderr" %in% output.param & length(output.param)==1){
    e <- new("ExpressionSet", assayData=assayDataNew(exprs=output$exprs, se.exprs=output$stderr), annotation=annotation(object))
  } else{
    if("stderr" %in% output.param){
      e <- new("frmaExpressionSet", assayData=assayDataNew(exprs=output$exprs, se.exprs=output$stderr), annotation=annotation(object))
    } else{
      e <- new("frmaExpressionSet", exprs=output$exprs, annotation=annotation(object))
    }
    if("weights" %in% output.param) e@weights <- output$weights
    if("residuals" %in% output.param) e@residuals <- output$residuals
  }
  
  return(e)
}
  

frmaBatch1 <- function(object, background, normalize, input.vecs, output.param, verbose){

  cdfname <- cleancdfname(cdfName(object))
  cdfname<-"hgu133plus2"
	platform <- gsub("cdf","",cdfname)
if(background == "rma"){
    if(verbose) message("Background Correcting ...\n")
    object <- bg.correct.rma(object)
  }

  if(is.null(input.vecs$normVec) | is.null(input.vecs$probeVec) | is.null(input.vecs$probeVarWithin) | is.null(input.vecs$probeVarBetween)){
    pkg <- paste(platform, "frmavecs", sep="")
    require(pkg, character.only=TRUE, quietly=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=eval(pkg))

    if(is.null(input.vecs$normVec)) input.vecs$normVec <- get(pkg)$normVec
    if(is.null(input.vecs$probeVec)) input.vecs$probeVec <- get(pkg)$probeVec
    if(is.null(input.vecs$probeVarWithin)) input.vecs$probeVarWithin <- get(pkg)$probeVarWithin
    if(is.null(input.vecs$probeVarBetween)) input.vecs$probeVarBetween <- get(pkg)$probeVarBetween
  }

  if(normalize == "quantile"){
    if(verbose) message("Normalizing ...\n")
    pm(object) <- normalize.quantiles.use.target(pm(object), input.vecs$normVec)
  }

  if(verbose) message("Summarizing ...\n")
  pns <- probeNames(object)
  pms <- log2(pm(object))
  
  tmp <- split(data.frame(pms, input.vecs$probeVec, input.vecs$probeVarWithin, input.vecs$probeVarBetween), pns)
  fit <- lapply(tmp, batchFit)
  
  exprs <- matrix(unlist(lapply(fit, function(x) x$exprs)), ncol=ncol(pms), byrow=TRUE)
  rownames(exprs) <- names(fit)
  colnames(exprs) <- colnames(pms)

  if("weights" %in% output.param){
    weights <- matrix(unlist(lapply(fit, function(x) x$w)), ncol=ncol(pms), byrow=TRUE)
    rownames(weights) <- pns
    colnames(weights) <- sampleNames(object)
  } else weights <- NULL

  if("stderr" %in% output.param){
    stderr <- unlist(lapply(fit, function(x) x$se))
    names(stderr) <- names(fit)
  } else stderr <- NULL

  if("residuals" %in% output.param){
    residuals <- apply(exprs,2,function(x) rep(x, table(pns)))
    residuals <- (pms-input.vecs$probeVec) - residuals
    rownames(residuals) <- pns
    colnames(residuals) <- sampleNames(object)
  } else residuals <- NULL
  
  colnames(exprs) <- sampleNames(object)
  
  return(list(exprs=exprs, stderr=stderr, weights=weights, residuals=residuals))
}

batchFit <- function(x){
  y.new <- as.matrix(x[,1:(ncol(x)-3)]) - x$input.vecs.probeVec
  yy.new <- as.vector(t(y.new))
  X <- matrix(rep(diag(ncol(y.new)),nrow(y.new)), nrow=length(yy.new), byrow=TRUE)
  Z <- diag(nrow(y.new))
  Z <- Z[rep(seq(nrow(Z)), each=ncol(y.new)),]
  G <- diag(x$input.vecs.probeVarBetween)
  R <- diag(rep(x$input.vecs.probeVarWithin, each=ncol(y.new)))
  V <- Z%*%G%*%t(Z) + R
  Vinv <- solve(V)
  e <- eigen(Vinv)
  C <- (e$vectors)%*%diag(sqrt(e$values))%*%t(e$vectors)
  yy.trans <- C%*%yy.new
  x.trans <- C%*%X

  fit.batch <- rlm(yy.trans ~ -1 + x.trans, maxit=100)
  b <- fit.batch$coef
  names(b) <- NULL
  u <- G%*%t(Z)%*%Vinv%*%(yy.new-X%*%b)
  offset.u <- mean(u)
  exprs <- b + offset.u
  w <- fit.batch$w
  se <- fit.batch$s
  return(list(exprs=exprs, w=w, se=se))
}
  


object1 <- frma1(data, summarize = "batch")

write.exprs(object1,"created_cdf_hgua_probe_gene.txt",sep="\t")
